Numerical analysis of impact processes of granular jets 
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Abstract. The rheology of a three-dimensional granular jet during an impact is investigated numerically. The cone-like 
scattering pattern and the sheet-like pattern observed in an experiment [X. Cheng, et al. Phys. Rev. Lett. 99, 188001 (2007)] can 
be reproduced through our calculation. We discuss the constitutive equation for granular jet impact in terms of our simulation. 
From the analysis of an effective friction constant, which is the ratio between the shear stress and the pressure the assumption 
of the zero yield stress would be natural in our setup and the shear visocity is not small in contrast to the suggestion by the 
experiment. 
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INTRODUCTION 

Impact processes of granular flow can be found in wide 
length scale, not only as problems of natural science but 
also as those of industrial applications [1, 2, 3, 4, 5, 6, 
7, 8]. A familiar application would be an ink-jet printing, 
which is an impact process of cohesive grains whose size 
are within nano-scale[2]. Recent experimental and theo- 
retical studies revealed interesting aspects of the impact 
processes of the granular flow. The impact of a granular 
flow onto a wall, produces a shock, which quantitatively 
agrees with the Mach cone produced by supersonic gas 
flow, at low volume fraction[3, 4, 5]. The impact dynam- 
ics of grains is also important as the geophysical prob- 
lems such as the formation of craters[6, 7, 8]. One of 
the recent interesting topics for the impact of the granu- 
lar flow would be the correspondence between granular 
flow and Quark Gluon Plasma(QGP), which is expected 
to behave like the fluid with very small shear viscosity, 
has been reported experimentally [1]. 

Recently, we report that the shear viscosity during 
the impact is not anomalous, although small shear stress 
observed in the experiment is reproduced through three- 
dimensional (3D) simulation. Thus, the correspondence 
between granular flow and QGP would be superficial [9]. 
However, through this analysis, zero yield stress, which 
is the residual stress without deformation, is assumed, 
though we presented three indirect evidences to support 
the assumption. In general, this assumption is a strong 
one, in particular, for the frictional case. 

In this paper, by introducing the effective friction coef- 
ficient and the inertia number, which are conventionally 
used for the dense granular flow[10, 11], the rheology 
of the granular jet is discussed. Here, we report that the 
shear viscosity consistent with our previous results, ex- 
cept its density dependence, and the assumption of the 



zero yield stress would be natural, by performing three- 
dimensional simulation [9]. 



MODEL 

We adopt the discrete element method (DEM) for mono- 
disperse soft core particles of the diameter rf[12]. When 
the particle i at the position r,- and the particle j at 
r,- are in contact, the normal force F-] is described as 
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the velocity v, of the particle i. The tangential force 
is given by F\- = min{|^ ; |,/xfj"}sgn(f|' / ), where ji is 
the local friction constant between contacting grains, 
sgn(x) = 1 for x > and sgn(x) = — 1 for otherwise, 
F-j = k, 8]j — rj, 8\j with the tangential overlap 5/ ; between 
i and j particles and the tangential component of relative 
velocity <5/ ; between i th and j th particles. Here, we 
adopt parameters k, — 0.2£„,rj, = 0.5r/„, ji = 0.2, k„ = 
4.98 x W 2 mul/d 2 , rj„ = 2.S8u /d and fi = 0.2, with 
incident velocity uq and the particle mass m. This set 
of parameters implies that the restitution coefficient for 
normal impact is e = 0.75 and duration time is t c = 
0.10d/uo. The value of ji is close to the experimental 
value for nylon spheres. [13] We adopt the second-order 
Adams-Bashforth method for the time integration with 
the time interval At = 0.02f c . 

Initial configurations are generated as follows: We pre- 
pare fee crystals and remove particles randomly to reach 
the desired density. We control the initial volume frac- 
tion 0o/0fcc = 00 before the impact as 0.30 < 0o < 0.90 
with volume fraction for a fee crystal 0f cc ~ 0.74 and 
20,000 particles are used. The initial granular temper- 
ature, which represents the fluctuation of particle mo- 



FIGURE 1. The side view of snapshots for the three dimensional simulation. For D tsI /Dj et = 0.8, cone-like scattered jet is 
reproduced (a), where black wall particles are hidden by scattered jet. For Aar/Djet = 2.5, sheet-like scattered jet is reproduced (b). 



tion, is zero. The wall consists of one-layer of parti- 
cles, which are connected to each other and with their 
own initial positions via the spring and the dashpot with 
spring constant k p = lO.Omu^/ d 2 and the dashpot con- 
stant r\ p = 5.077„, respectively. 

Experimentally, it is known that the scattered state ex- 
hibits the crossover from a cone-like pattern and a sheet- 
like pattern by changing Aar/Aet with the jet diameter 
£>j et and the target diameter D tar . The crossover can be 
reproduced through DEM, where the jet diameter Dj et is 
fixed Dj et /d — 4.5 (Fig. 1). White particles (open circles 
in Fig.l) are grains and black solid ones are wall parti- 
cles. The Figure 1 (a) is a typical con-like pattern with 
Aar/Aet = 0-8. The Figure 1 (b) is an example of the 
sheet-like pattern with Aar/Aet = 2-5 (b). We note that 
wall particles are hidden in (a). 

RHEOLOGY OF GRANULAR JETS 

We evaluate physical quantities near the wall at the 
height z = Az — 5.0d from the wall z — 0. D- jet /d — 
10.0 and D tar /d = 22.0 . We adopt the cylidrical co- 
ordinate whose symmetric axis is chosen to be the jet 
axis, and divide calculation region into the radial direc- 
tion r = 0,Ar, ■ ■ ■ ,5Ar, with Ar = /?tar/5 with the target 
radius R tar . Then we estimate physical quantities in the 
corresponding mesh region with kAr < r < (k+1 )Ar 
(k = 0, 1, • • • ,5), where r is the distance from the sym- 
metric axis of the cylindrical coordinate. 

We calculate stress tensor as in Ref. [14]. The micro- 
scopic definition of the stress tensor at r is given by 

(r ) = \ £ muiaUjp + ^ £ F% r'J , ( 1 ) 



where i and j are indices of particles, a,j3 = r, d,z de- 
notes cylindrical coordinates and denotes the summa- 
tion over the particles denoted by ;' located at r. Here, z 
axis is parallel to the incident jet axis, and V is the vol- 
ume of each mesh at r and w (Ce (r) = v' a — v a (r) with the 
mean velocity v«(r) in the mesh at r. To calculate the 
stress tensor in cylindrical coordinates, we firstly calcu- 
late O a ipi in Cartesian coordinate, a', J3' =x,y,z, whose 
origin is the same as cylindrical one, and transform it into 
that for cylindrical one. 



Velocity profile 

The profile of v r and y/T g /m is shown in Fig. 2, 
with the granular temperature T g (r) = £ ice mM? a (r)/3Af. 
Ellowitz et al. suggests that the dead zone, where the 
motion of grains is frozen, exist near the target in two 
dimension (2D) [15, 16]. However, as is shown in Fig. 
2, although the velocity of grains at the center is small, 
the fluctuation of the particle velocity y/T g /m is the 
largest at the center. Thus, the motion of particles near 
the target in 3D is not frozen. Namely, there is no dead 
zone in 3D granular jets. It should be noted that, in our 
2D calculation, T g is small at the center, i.e. the dead zone 
actually exists, which will be reported elsewhere. 



Pressure 

Following hydrodynamical model for the granular 
flow proposed by Garzo and Dufty[17], pressure P = 
La a aa /3 is conventionally given by 

-L = l+2<l>(l+e) X , (2) 
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FIGURE 2. The profile of the radial component of the veloc- 
ity field v r and ^T g /m. 
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FIGURE 3. Analysis on the /i* vs / plane. Numerical datum 
can be fitted into black solid lines /i* = al. 
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where 0/ = 0.49, f = 0.64, number density n, vol- 
ume fraction = nncP /6 and radial distribution function 
*[18]. 

For the frictional case, in general, five equations for 
rotational degree of freedom are necessary, in addition to 
those for the translational one. However, ten equations 
for frictional grains can be reduced to five equations 
by introducing effective restitution coefficient e, if the 
friction constant fi is small [19, 20, 21]. According to this 
simplification we use the effective restitution coefficient 
e = 0.616 for e = 0.75 and jj, = 0.2, for frictional case in 
the following analysis. 

In our setup, the empirical relation (2) gives a good 
approximation for < <j>f, while the deviation between 
numerical data and theoretical curve exists for denser re- 
gions near the symmetric axis which may result from the 
singularity r ~ of the cylindrical coordinate. For later 
analysis, thus, we adopt the equation of state P((j>,T g ) = 
nT g {l+2(j)(l+e)x}[9]. 



Friction coefficient 

Let us analyze the effective friction coefficient for 
macroscopic motion of a collection of grains, fol- 
lowing the ref. [22]. We estimate strain rate D n as 

dv r (r,Az/2)/dz £ii (v r (r,3Az/4) - v r (r,Az/4))/(Az/2) 
and dv z (r,z)/dr ~ (v z (r+Ar/2,z) -v z (r- Ar/2,z))/Ar. 
Since we evaluate the physical quantities near 
the wall, the mesh < z < Az is divided into 
< z < Az/2 and Az/2 < z < Az to calculate 
dv r (r,Az/2)/dz and < r < Rt m - is divided into 
< r < Ar/2, Ar/2 < r < 3Ar/2, • • •. 



Introducing the effective friction co efficien t jj.* = 
— <7 n /P and the inertia number / = D rz ^JP /md, we plot 
the observed data fi* vs / for several 0o in Fig. 3. From 
the result of our simulation, the obtained effective fric- 
tion coefficient ji* can be fitted by ji* — al with a con- 
stant a within error bars, where fitting values are a — 
0.240 and a = 0.223 for frictional and frictionless case, 
respectively. The solid lines represent the corresponding 
fitting lines. Judging from the fitting, the assumption of 
zero yield stress in our setup[9] would be natural. 

Through this analysis, the constitutive equations are 
obtained as 
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1* = — V60(l+20(l+ e )z). (5) 

Here, non-dimensional shear viscosity tj* = 77 * (0 ) is 
introduced as tj* = TJ /tjo with the shear viscosity tj and 
Tjo = 5y / 7tmT g /16d 2 by convention [17]. 

The obtained shear viscosity is proportional to v/5p, 
which is consistent with our previous results [9]. How- 
ever, the obtained density dependence of tj* is differ- 
ent from ref. [17]. The non-dimensional shear viscosity 
based on the kinetic theory Tjj! in (0) is compared with our 
result in Fig. 4. The solid line and the dashed line de- 
notes for the frictional case and the frictionless case, re- 
spectively. The obtained shear viscosity is less than 80% 
of that of the kinetic theory. It should be noted that the 
obtained shear viscosity is also finite, even through the 
analysis based on the effective friction coefficient. 



DISCUSSION 

The difference of density dependence of non- 
dimensional shear viscosity would be understood in 
the followings. In the analysis on pL* vs / plane, the 
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FIGURE 4. Comparison of the non-dimensional shear vis- 
cosity based on the kinetic theory '?y n (<^)[17] with our result, 
where the solid line and the dashed line denotes for the fric- 
tional case and the frictionless case, respectively. 



density dependence of transport coefficient is assumed 
to appear only through pressure P, i.e. 

7](P,T g ,<j>) = T 1 (P(4>,T g )). (6) 

Thus, once the equation of state is determined, 77 * (0 ) is 
uniquely obtained. On the other hand, in ref. [17], be- 
cause density dependence of 77j* in results from the ki- 
netic theory, results are different. At present, we could 
not judge which viscosity is better. However, we should 
note that the tj °< ^JT^ can be obtained, even under the 
assumption of eq. (6). 

We comment here that fi* — in two dimensions 
cannot be fitted by fi* — al. Because grains are easily 
packed through the impact in 2D, dense flow and related 
frictional phenomena, where fi* would be greater by 
a degree of magnitude than that in 3D, emerges. The 
existence of the dead zone may cause the difference 
between 2D and 3D results. The frictional phenomena 
in 2D will be reported elsewhere. 

CONCLUSION 

We have numerically investigated the constitutive equa- 
tion for the granular jet impact in 3D, introducing the 
analysis on fi* vs / plane. The dead zone, which is sug- 
gested by Ellowitz et a/[16], can not be reproduced in 
our 3D study, although the velocity field at the center is 
small. Rheological results are consistent with our previ- 
ous results [9], except for the density dependence of the 
shear viscosity, which results from the assumption in Eq. 
(6). Judging from the analysis, the assumption of zero 
yield stress would be natural[9]. 
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